Última actualización: 25 mayo 2020
En esta entrega se describe la estrategia de calibración usada con 5 estaciones de caudal. Los datos observados de caudal han sido previamente tratados (ver: Tratamiento de caudales).
La siguiente tabla muestra los periodos de validación.
t23 = c("2004-01-01","2005-12-31")
t30 = c("2005-01-01","2006-12-31")
t36 = c("2008-01-01","2009-12-31")
t38 = c("2007-01-01","2008-12-31")
t41 = c("2006-01-01","2007-12-31")
T.valida = rbind(t23,t30,t36,t38,t41)
kable(T.valida)| t23 | 2004-01-01 | 2005-12-31 |
| t30 | 2005-01-01 | 2006-12-31 |
| t36 | 2008-01-01 | 2009-12-31 |
| t38 | 2007-01-01 | 2008-12-31 |
| t41 | 2006-01-01 | 2007-12-31 |
Las siguientes figuras muestran en gris los periodos de validación, el resto de la serie se utiliza para calibración.
#lectura de datos
qobs <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/qobs_5est_2004_2012.RDS")
# ajuste de caudal
qobs$q23 = 1.630*qobs$q23
qobs$q30 = 0.789*qobs$q30
qobs$q36 = 1.184*qobs$q36
qobs$q38 = 0.724*qobs$q38
qobs$q41 = 1.120*qobs$q41
#conversion a mm
qobs$q23 = convertFlow(qobs$q23, from="cumecs", to="mm", area.km2 = 687)
qobs$q30 = convertFlow(qobs$q30, from="cumecs", to="mm", area.km2 = 1092)
qobs$q36 = convertFlow(qobs$q36, from="cumecs", to="mm", area.km2 = 2744)
qobs$q38 = convertFlow(qobs$q38, from="cumecs", to="mm", area.km2 = 3159)
qobs$q41 = convertFlow(qobs$q41, from="cumecs", to="mm", area.km2 = 4896)fig = plotly_empty()
fig = fig %>% add_segments(x = ~t23[1], y=max(qobs$q23,na.rm = T), yend = max(qobs$q23,na.rm = T),
xend = t23[2], name="Validación Paso Pache",
line = list(color = 'black', width = 2))
fig = fig %>% add_segments(x = ~t30[1], y=max(qobs$q30,na.rm = T), yend = max(qobs$q30,na.rm = T),
xend = t30[2], name="Validación Paso Roldán",
line = list(color = 'blue', width = 2))
fig = fig %>% add_segments(x = ~t36[1], y=max(qobs$q36,na.rm = T), yend = max(qobs$q36,na.rm = T),
xend = t36[2], name="Validación Fray Marcos",
line = list(color = 'red', width = 2))
fig = fig %>% add_segments(x = ~t38[1], y=max(qobs$q38,na.rm = T), yend = max(qobs$q38,na.rm = T),
xend = t38[2], name="Validación San Ramon",
line = list(color = 'orange', width = 2))
fig = fig %>% add_segments(x = ~t41[1], y=max(qobs$q41,na.rm = T), yend = max(qobs$q41,na.rm = T),
xend = t41[2], name="Validación Paso Pache",
line = list(color = 'green', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q23, type="scatter",mode = "lines",
name="Paso de los Troncos", line = list(color = 'black', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q30, type="scatter",mode = "lines",
name="Paso Roldán", line = list(color = 'blue', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q36, type="scatter",mode = "lines",
name="Fray Marcos", line = list(color = 'red', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q38, type="scatter",mode = "lines",
name="Paso San Ramón", line = list(color = 'orange', width = 2))
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q41, type="scatter",mode = "lines",
name="Paso Pache", line = list(color = 'green', width = 2))
fig <- fig %>% layout(#legend = list(x = 1, y = 0.9),
yaxis = list(title="Runoff (mm)",
showticklabels = TRUE, ticks="outside",
showline=TRUE, showgrid=TRUE),
xaxis = list(title="",showticklabels = TRUE,
ticks="outside",showline=TRUE,
showgrid=TRUE, range=as.Date(c("2004-01-01", "2012-12-31"))),
height=300,
width=700)
figrm(list = ls())
# Librerias
library(SWATplusR)
library(stringr)
library(hydroGOF)
library(hydromad)
library(tibble)
library(hydroPSO)
# Caudales observados
obs.data <- readRDS("./data/qobs_5est_2004_2012_para_calibracion.RDS")
# ajuste de caudal
obs.data$q23 = 1.630*obs.data$q23
obs.data$q30 = 0.789*obs.data$q30
obs.data$q36 = 1.184*obs.data$q36
obs.data$q38 = 0.724*obs.data$q38
obs.data$q41 = 1.120*obs.data$q41
# SWAT project mensual
swat_project = "./SWAT_SL_20200401/"
# Funcion de optimización
swat_optim <- function(mypar, obs.data, project_dir, start_sim, end_sim,runrun_path, head_log) {
names(mypar) = head_log[25:length(head_log)]
path_run = paste0(runrun_path,"run_",Sys.getpid())
par_set = tibble("CN2::CN2.mgt|change = pctchg" = mypar["CN2"]*100,
"ESCO::ESCO.hru|change = absval" = mypar["ESCO"],
"SOL_AWC::SOL_AWC.sol|change = pctchg " = mypar["SOL_AWC"]*100,
"ALPHA_BF::ALPHA_BF.gw|change = absval" = mypar["ALPHA_BF"],
"GWQMN::GWQMN.gw|change = absval" = mypar["GWQMN"],
"GW_DELAY::GW_DELAY.gw|change = absval |" = mypar["GW_DELAY"],
"REVAPMN::REVAPMN.gw|change = absval" = mypar["REVAPMN"],
"GW_REVAP::GW_REVAP.gw|change = absval" = mypar["GW_REVAP"],
"OV_N::OV_N.hru|change = absval" = mypar["OV_N"],
"SLSUBBSN::SLSUBBSN.hru|change = pctchg" = mypar["SLSUBBSN"]*100,
"HRU_SLP::HRU_SLP.hru|change= pctchg" = mypar["HRU_SLP"]*100,
"USLE_K::USLE_K.sol|change= pctchg" = mypar["USLE_K"]*100,
"USLE_P::USLE_P.mgt|change= pctchg"=mypar["USLE_P"]*100,
"P_UPDIS::P_UPDIS.bsn|change= absval"= mypar["P_UPDIS"],
"FRT_SURFACE::FRT_SURFACE.mgt|change=absval" = mypar["FRT_SURFACE"],
"FRT_KG::FRT_KG.mgt|change=pctchg" = mypar["FRT_KG"]*100,
"NPERCO::NPERCO.bsn|change= absval"=mypar["NPERCO"],
"ERORGN::ERORGN.hru|change= absval"= mypar["ERORGN"],
"PPERCO::PPERCO.bsn|change= absval"= mypar["PPERCO"],
"CDN::CDN.bsn|change= absval"= mypar["CDN"],
"CMN::CMN.bsn|change= absval"= mypar["CMN"],
"PSP::PSP.bsn|change= absval"= mypar["PSP"]
)
# run swat
qsim = run_swat2012(project_path = project_dir,
output = list(FLOW = define_output(file = "rch",
variable = "FLOW_OUT",
unit = c(23,30,36,38,41)),
PT = define_output(file = "rch",
variable = "TOT_Pkg",
unit = c(23,41)),
NO3 = define_output(file = "rch",
variable = "NO3_OUT",
unit = c(23,41))),
parameter = par_set,
start_date = as.Date("2000-01-01"),
end_date = as.Date("2012-12-31"),
years_skip = 4,
run_path=path_run,
quiet = TRUE,
output_interval="d",
add_parameter = FALSE)
# objective function
nse23 = NSE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
nse30 = NSE(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
nse36 = NSE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
nse38 = NSE(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
nse41 = NSE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
kge23 = KGE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
kge30 = KGE(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
kge36 = KGE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
kge38 = KGE(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
kge41 = KGE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
pbias23 = pbias(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
pbias30 = pbias(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
pbias36 = pbias(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
pbias38 = pbias(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
pbias41 = pbias(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
# objective quality
qlt_41 = data.frame(flow_lps_km2=1000*qsim$FLOW_41/4896,
PT_mgs_km2 = 1000000*qsim$PT_41/(24*60*60*4896),
NO3_mgs_km2 = 1000000*qsim$NO3_41/(24*60*60*4896),
qobs_lps_km2 = 1000*obs.data$q41/4896,
PTobs_mgs_km2 = 0.16*(1000*obs.data$q41/4896)^1.02,
NO3obs_mgs_km2 = 0.1*(1000*obs.data$q41/4896)^1.41)
qlt_41[qlt_41<=0] = NA
qlt_41 = na.omit(qlt_41)
qlt_23 = data.frame(flow_lps_km2=1000*qsim$FLOW_23/687,
PT_mgs_km2 = 1000000*qsim$PT_23/(24*60*60*687),
NO3_mgs_km2 = 1000000*qsim$NO3_23/(24*60*60*687),
qobs_lps_km2 = 1000*obs.data$q23/687,
PTobs_mgs_km2 = 0.06*(1000*obs.data$q23/687)^0.93,
NO3obs_mgs_km2 = 0.04*(1000*obs.data$q23/687)^1.41)
qlt_23[qlt_23<=0] = NA
qlt_23 = na.omit(qlt_23)
fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_41$PT_fit = a*(qlt_41$flow_lps_km2)^b
fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_41$NO3_fit = a*(qlt_41$flow_lps_km2)^b
fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_23$PT_fit = a*(qlt_23$flow_lps_km2)^b
fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]
qlt_23$NO3_fit = a*(qlt_23$flow_lps_km2)^b
kge_PT_41 = KGE(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
kge_NO3_41 = KGE(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
kge_PT_23 = KGE(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
kge_NO3_23 = KGE(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
pbias_PT_23 = pbias(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
pbias_PT_41 = pbias(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
pbias_NO3_23 = pbias(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
pbias_NO3_41 = pbias(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
FO = 0.5*(kge23 + kge30 + kge36 + kge38 + kge41)/5 + 0.25*(kge_PT_41 + kge_PT_23)/2 + 0.25*(kge_NO3_41 + kge_NO3_23)/2
# logfile
vout = data.frame(FO,kge23,kge30,kge36,kge38,kge41,
nse23,nse30,nse36,nse38,nse41,
pbias23,pbias30,pbias36,pbias38,pbias41,
kge_PT_23, kge_PT_41, kge_NO3_41, kge_NO3_23,
pbias_PT_23, pbias_PT_41, pbias_NO3_23, pbias_NO3_41,
t(mypar))
colnames(vout) = head_log
log.file = paste0(runrun_path,"log_",Sys.getpid(),".txt")
if(!any(list.files(runrun_path)==paste0("log_",Sys.getpid(),".txt"))){
write.table(t(head_log), log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
}
ff = read.table(log.file, col.names = head_log,colClasses = "character")
ff = rbind(ff,vout)
write.table(ff, log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
#output
return(FO)
}
CN2 = c(-0.2,0.2)
ESCO = c(0.1,0.45)
SOL_AWC = c(-0.2,.2)
ALPHA_BF = c(0.4,0.8)
GWQMN = c(800,4000)
GW_DELAY = c(20,45)
REVAPMN = c(300,700)
GW_REVAP = c(0.02, 0.01)
OV_N = c(0.4, 0.6)
SLSUBBSN = c(-0.1, 0.3)
HRU_SLP = c(-0.1,0.1)
USLE_K = c(-0.2,0.2)
USLE_P = c(-0.2,0.2)
P_UPDIS = c(20, 90)
FRT_SURFACE = c(0.1,0.6)
FRT_KG = c(-0.3,0.3)
NPERCO = c(0.5,1)
ERORGN = c(2,4)
PPERCO = c(10,17.5) #bsn absval
CDN = c(0,3) #bsn absval
CMN = c(0,0.03) #bsn absval
PSP = c(0.01,0.8)
names_par = c("CN2","ESCO","SOL_AWC",
"ALPHA_BF","GWQMN","GW_DELAY","REVAPMN","GW_REVAP","OV_N","SLSUBBSN","HRU_SLP",
"USLE_K","USLE_P","P_UPDIS","FRT_SURFACE","FRT_KG","NPERCO","ERORGN","PPERCO","CDN","CMN","PSP")
# CN2 ESCO SOL_AWC ALPHA_BF GWQMN GW_DELAY REVAPMN GW_REVAP OV_N SLSUBBSN
par_lwr = c(CN2[1],
ESCO[1], SOL_AWC[1], ALPHA_BF[1], GWQMN[1], GW_DELAY[1], REVAPMN[1], GW_REVAP[1], OV_N[1],SLSUBBSN[1],
HRU_SLP[1], USLE_K[1], USLE_P[1], P_UPDIS[1], FRT_SURFACE[1], FRT_KG[1], NPERCO[1], ERORGN[1], PPERCO[1], CDN[1], CMN[1],PSP[1])
par_upr = c(CN2[2],
ESCO[2], SOL_AWC[2], ALPHA_BF[2], GWQMN[2], GW_DELAY[2], REVAPMN[2], GW_REVAP[2], OV_N[2],SLSUBBSN[2],
HRU_SLP[2], USLE_K[2], USLE_P[2], P_UPDIS[2], FRT_SURFACE[2], FRT_KG[2], NPERCO[2], ERORGN[2], PPERCO[2], CDN[2], CMN[2],PSP[2])
par_init = (par_upr + par_lwr)/2
# Optimización
head_log = c("FO","KGE23","KGE30","KGE36", "KGE38","KGE41",
"NSE23","NSE30","NSE36","NSE38", "NSE41",
"pbias23","pbias30","pbias36","pbias38", "pbias41",
"kge_PT_23", "kge_PT_41", "kge_NO3_41", "kge_NO3_23",
"pbias_PT_23", "pbias_PT_41", "pbias_NO3_23", "pbias_NO3_41",
names_par)
fr = list.files("./run_optim_mensual/", pattern = "*.txt", full.names=TRUE)
file.remove(fr)
opt_pso = hydroPSO(par=par_init, fn = swat_optim,lower=par_lwr, upper=par_upr,
control = list(MinMax = "max",
maxit=500,
reltol=0.00001,
parallel="parallelWin",
par.nnodes=12,
par.pkgs = list("hydromad",
"SWATplusR",
"stringr",
"hydroGOF",
"tibble"),
write2disk=FALSE,
REPORT=1,
npart=48,
normalise=TRUE),
obs.data = obs.data,
project_dir = swat_project,
runrun_path = "./run_optim_mensual/",
head_log = head_log
)A work by Rafael Navas